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The method of a determination of the Primary Cosmic Ray mass composition is presented. Data 
processing is based on the theoretical model representing the integral muon multiplicity spectrum 
as the superposition of the spectra corresponding to different kinds of primary nuclei. The method 
consists of two stages. At the first stage, the permissible intervals of primary nuclei fractions fi 
are determined on the base of the EAS spectrum vs the total number of muons {E^ > 235 GeV). 
At the second stage, the permissible intervals of fi are narrowed by fitting procedure. We use the 
experimental data on high multiplicity muon events (n^j > 114) collected at the Baksan underground 
scintillation telescope. Within the framework of three components (protons, helium and heavy 
nuclei), the mass composition in the region 10^^ — 10^® eV has been defined: fp — 0.235 ± 0.02, 
fHe = 0.290 ± 0.02, fn = 0.475 ± 0.03. 

PACS numbers; 26.40.+r 



I. INTRODUCTION 

Up to now the energy spectrum and the mass composition of primary cosmic rays (CR) have been measured by 
direct methods (on satellites or stratosphere balloons) up to energies of about ~ 100 TeV (i?„ is the primary 
nucleus energy). At higher energies the information about CR is obtained with the help of indirect methods which 
consist in a measurement of different parameters of extensive air showers (EAS). Parameters relating to the energy 
spectrum are in rather good agreement between each other [1 — 11], while the data on the mass composition are 
inconsistent enough (see table H]). 

Table I: Some parameters of the CR energy spectrum and mass composition obtained in different experiments. The second 
column shows the publish year and reference. 



Detector 


Reference iJfe,PeV 


71 


72 


< InA >,{p + a) 


MSU 


97, fl, 2] 


3 


2.7± 


3.1± 


0.62 


0.24 {p + a) 


EAS - TOP 


98, [3] 


3 


2.76 ± 0.03 


3.19 ±0.06 


7.1 - 


8.9 {A,ff) 


HEGRA 


97, [4] 


2 


2.60 ±0.10 


3.00 ±0.15 


0.60 


^ 0.45 {p + a) 


HEGRA 


98, f6] 




2.60 ±0.10 


3.00 ±0.15 


2.0 - 


^ 2.4 (< InA >) 


HEGRA 


99, f7] 


3.4 


2.67 ±0.03 


3.33 ± 0.40 


0.55 


0.48 (p ± a) 


KASCADE 


99, [8] 


4 


2.7 


3.1 


0.70 


0.50 (p + a) 


CASA - MIA 


99, [9] 


3 


2.66 ±0.02 


3.00 ±0.05 


1.3 - 


3 (< InA >) 


DICE 


00,flO] 


3 


2.66 ±0.02 


3.00 ±0.05 


1.5 - 


0.9 (< InA >) 












0.60 


- 0.50 - 0.65 (p + a) 


CASA - 


01, fU] 


2-3 


2.72 ±0.02 


2.95 ±0.02 


0.62 


- 0.80 - 0.43 (p ± a) 


- BLANCA 










1.7 - 


1.0 - 1.4 (< InA >) 



To investigate the CR mass composition, the EAS parameters are measured which have to be distinct for EAS 
produced by different kinds of nuclei. These are the number of muons in EAS , the maximum depth in atmosphere 
Xm, fluctuations of the maximum depth a{Xm), a steepness of lateral distribution of particles near EAS core p(r) 
etc. 

The main problem of indirect methods of CR investigation is that the information about both the energy spectrum 
and the mass composition must be obtained from the same data sample. This leads hereto that the determination 
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of CR mass composition is very difficult problem. There are numerous and various uncertainties related with an 
interaction model and methods of measurement of EAS characteristics (having, as a rule, considerable errors) and 
with the CR energy spectrum. These difhculties lead to large spread of obtained results. In Table HI the parameters 
of the CR energy spectrum (Ek is the "knee" energy, 71 and 72 are slope exponents of the spectrum before and after 
the "knee") and the mass composition (< InA > is the average logarithm of the number of nucleons in a nucleus, 
(p + a) is the light nuclei fraction) obtained in different experiments are presented. One can see that the trend to 
weighting mass composition is violated by results reported in [lol [Tl| . 

In the paper we present a method of the determination of CR mass composition on the base of data on the EAS 
spectrum vs the total number of high energy muons J(n^). 

The paper is constructed as follows. In Section 2, entry conditions are described. In Section 3, we present the 
method of the determination of primary nuclei fractions intervals A/^, which ensure an agreement with experimental 
data within the limits of one standard deviation. The realization of the method is shown in Sections 4 and 5. In 
Section 6, intervals of Afi determined at the first stage are narrowed by fitting procedure. Sections 7 and 8 are 
Discussion and Conclusion. 



II. INITIAL CONDITIONS 



We shall use the data on high multiplicity muon events (n^ > 114) collected at the Baksan underground scintillation 



telescope 12|. In [13| the muon multiplicity spectrum (i.e., the number m of muons hitting the facility at unknown 
position of EAS axis) at m > 20 was measured at zenith angles 6 < 20°. The threshold energy of muons coming from 
this solid angle is 235 GeV. 

It is known, the muon multiplicity spectrum depends on the facility geometry and selection conditions of the 
experiment. This leads to that multiplicity spectra cannot be compared with each other. It would be better to 
present the data as a function of some invariant variable which does not depend on experiment conditions. In our 
opinion, the total number of muons in the EAS, n^, can be chosen as a such variable. 

In papers [T3 - [T6| we developed the method of recalculation from multiplicity spectrum to the EAS spectrum vs 
the total number of muons, I{n^). The formulation of the task is following: let F(m) be the integral multiplicity 
spectrum obtained at a certain facility. Let us define the parameter A(to) — mi/n^, which is the average fraction of 
muons hitting the facility in the case when the latter is crossed by mi > m muons. Assuming then = to/A(to), 
we will obtain the integral spectrum of EAS vs the total number of muons 

/(> n^) = -pJ-^F{^), (1) 
G{m) 

here G{m) is the acceptance of the facility for a collection of events with muon multiplicity > m. (It should be 
explained that m has to be high enough, for example m > 20.) 

The numerical values of parameters A{m) and G{m) calculated in [isl. [Tg} with regard to the real structure of the 
facility are presented in Table HIl N{> m) is the experimental number of events with the muon multiplicity > m [l3j . 
Nonmonotony of the N{> m) is due to the different exposure time. Tree- 



Table II: Numerical values of A(m) and G{m). The error in calculation of A(m) and G{m) does not exceed 2-3%. Non-integer 
values of m are obtained because of corrections at trajectories reconstruction. 



m 


N{> m) 


Tfec^ 10 s 


A(m) 


G(m), ■ sr 


rifj. 


21.9 


127 


1.67 


0.280 


60.5 


78.2 


32.9 


547 


19.33 


0.289 


57.9 


113.9 


44.5 


270 


19.33 


0.295 


56.6 


150.8 


56.5 


164 


19.33 


0.299 


54.8 


188.6 


82.1 


66 


19.33 


0.306 


53.2 


268.2 


124.9 


49 


41.93 


0.313 


51.6 


399.3 


211.6 


7 


41.93 


0.319 


50.4 


663.8 



The method developed is universal and allows to combine results obtained in different experiments with muon 
bundles. In our case we have combined the results reported in [l^ and [H, [l3| and obtained the EAS spectrum 
vs the total number of muons in the range 75 < < 4000, which corresponds to the primary energy range of 
10^^ < Ejs, < 10^^ eV (Fig.l). It should be clarified that the data at > 2000 are obtained for the muon threshold 
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Figure 1: Squares are the EAS spectrum vs riu (experimental data). The muon threshold energy is Eth = 235 GeV if < 
1000 and Eth ~ 220 GeV at > 1000 [l4l. [la]. Crosses show the muon multiplicity spectrum obtained in {m and F{m) 
correspond to the multiplicity spectrum). Solid curves are expected fluxes {Eth ~ 220 GeV) for the case Ek^Z -3- W^'^ eV, 
dashed curves - the case Ek = S ■ 10^^ eV/nucleus. Dotted curves show expected fluxes for the case Ek = Z ■ 3 ■ 10^^ eV at 
Eth = 235 GeV. Numbers near curves denote the mass composition variants: 1 is the "standard" (low energy) composition, 
2 is the composition |[2]). 



energy Et}^ = 220 GeV, while the points at < 700 have Etfi — 235 GeV which is the threshold energy in the 
experiment [l^l- (We do not recalculate the data to the same threshold energy to avoid additional errors.) In Fig.l, 
the expected fluxes are calculated for Eth — 235 GeV (n^ < 1000, dotted curves) and Eth = 220 GeV (n^ > 1000, 
solid and dashed curves). Numbers near curves denote the mass composition variants: 1 is the low energy composition 
(the nuclei fractions in percentage are 39, 24, 13, 13, 11), 2 is the composition Note there is no normalization in 
Fig.l. 

The data at < 700 can be used for retrieval of information on the CR mass composition in the region E^^ = 
10^^ — 10^® eV. Let us remark here that the data at m = 124.9 and m = 211.6 in [3] were obtained with essential 
systematic errors: according to our estimates, the values of m in these points are underestimated 4% and 10% 
respectively jl^, therefore we restrict ourselves to the data at m < 82 (n^ < 270). 

As initial conditions, we use the mass composition obtained by Swordy [l^ with the help of compilation of results 
of direct measurements at energies ~ 100 TeV per nucleus^ (A is the number of nucleons in a nucleus) 

p He CNO Ne-S Fe 
A 1 4 14 28 56 (2) 
/,% 25 31 19 12 13 

and the proton flux at the energy Ep — 100 TeV measured in the JACEE experiment poj. 

Dp{100 TeV) = 2.95 x 10"^° {m^ ■ s ■ sr ■ GeVy^ (3) 



^ in comparison with the composition presented in [Tgl ]. in ((J} the proton fraction is increased by 5% (at the expense of helium nuclei) in 
accordance with data of |20|I 
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Then the total flux of nuclei with energy of £; = 100 TeV is equal to Aai = -Dp/0.25 = 11.8 xlQ-^" {m'^ ■ s ■ sr ■ GeV)-'^ , 
that is in a good agreement with the result obtained at Tibet array [2l|. In this case, the mass composition at the 
same energy per nucleon is 



F(100 TeV) 



( 0.88633 \ 
' 0.10412 ' 
0.007584 
0.001474 
V 0.000492 / 



(4) 



and the total flux of nuclei with the energy 100 TeV/nucleon is equal to 

13^(100 TeV) = ^DtotfjAj^ '^ = 3.329 x 10"^° {m^ ■ s ■ sr ■ GeV)-\ (5) 
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Our goal is to determine the mass composition evolution (from the composition ([2])) into the region — 10^^ — 10^^ eV 
on the base of data on the multiplicity of high energy muons (£'^ > 235 GeV) in EAS (see table|lTl). To this end, we will 
use the measured fluxes of multiple muon events with the multiplicity into differential intervals n^i < < 
At the first stage, we determine the permissible intervals of primary nuclei fractions fi which ensure an agreement 
with experimental data within the limits of one standard deviation (Sections 4 and 5). And then, we refine the results 
with the help of fitting procedure (Section 6). 

To obtain the more certain results we fix the CR energy spectrum, namely we adopt the conservative scenario: 

i) the slope change of the spectrum occurs at the same energy per unit charge Ek{Z) — 3 PeVxZ, 

ii) the spectra of all nuclei kinds have the slope exponents 71 — 2.7 before the "knee" and 72 = 3.1 after the "knee" 

Da{E) = IaE-'-'{1 + E/Ek{Z))-'-\ (6) 

As is seen from table |T] this scenario is supported by experimental data well enough. 

It should be emphasized, we do not attempt to use the data at > 2000 because the energy spectra of primary 
nuclei at £'„ > 10^^ eV are poorly understood. 

III. EQUATIONS 

The fiux of events with muon multiplicity n > produced by nuclei with A nucleons can be written in the form 

00 

lAin > rv) = J Da{E)Pa{E, n > n^)dE , (7) 

Eth(A) 

here Da{E) is the differential flux of nuclei of kind A, Pa{E, n > n^) is the probability that the number of muons 
(with Ef^ > 235 GeV) in EAS produced by nucleus "A" (with energy E per nucleon) is n > n^, Eth{A) is the threshold 
energy of nuclei with A nucleons. 

We assume that the multiplicity of muons in EAS is described by the negative binomial distribution Ba{E, n) 
(Appendix B), then 

Pa{E, n> n^) = Ba{E, n) 

n>n^ 

Taking into account that Da{E) — Dn{E) x Fa(£'), we rewrite ([7]) so 

00 

lAin > n^^) = J Fa{E)Dn{E)Pa{E, n > n^,)dE, (8) 

EthiA) 

and the flux of events with < ri^ < has the form 

JAi^rifj^i) = lA{n>n^^i) - lA{n>n^^i^i+i)) ^ J Fa{E)Dn{E)Pa{E, n> n^^i)dE- 

'-'"^ _ (9) 

- / Fa{E)Dn{E)Pa{E, n>n^,+l)dE = FAR^A, 

El+\A) 
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Figure 2: Energy distributions of protons and iron nuclei making a contribution to the flux of muon events with 114 < < 151. 
Areas under curves (p and Fe ) are equal to 1. The widths of distributions at half-height and fractions of events in these regions 
are indicated. 

where the first index of the matrix i?y points out to muon multiplicity (n^i) and the second one pertains to a nucleus 
sort. Fji is the fraction of nuclei "A" averaged over the energy region which gives the main contribution in the integral 
(as is seen in Figj^l the region is rather narrow). Thus we work in the approximation Fj(E) ~ Fj — const and 
will drop the symbol of averaging hereinafter. 

A fast decrease of the CR flux with the energy is an important simplifying factor. In consequence, the main 
contribution to the muon bundles flux at any threshold multiplicity is originated from nuclei whose energies are 
in a rather narrow region. In FiglJl the energy distributions of protons and iron nuclei making a contribution to the 
flux of events (EAS) with 114 < ^^(E'^ > 235 GeV) < 151 are shown. The widths of distributions at half-height are 
1745 - 1150 TeV for iron nuclei and 4985 - 2485 TeV for protons. 

To avoid possible methodical errors, we use only 4 points in the spectrum Itot{^ n^^) = /(> rt^): at ~ 114, 151, 
189, 268 (the point at = 78 has a different exposure time and we do not use it in the present work (see table |TT|). 
In table Hin the input data are presented: muon multiplicity intervals An^, the numbers and fluxes of events in given 
intervals of . 

Table III: Integral (7) and finite-difference (J) fluxes of events (EAS) with the given number of muons (_E^ > 235 GeV). The 
flux J(An,j) is defined according to ((9]). 





7(> Uf,) X 10^ (m^ ■ s • sr)-^ 




iV(An^) 


Jx 10^ (m^-s-sr)-^ 


114 


4.887 ± 0.209 


114 - 151 


277 


2.419 ±0.145 


151 


2.468 ±0.150 


151 - 189 


106 


0.920 ± 0.089 


189 


1.548 ±0.121 


189 - 268 


98 


0.906 ± 0.092 


268 


0.642 ± 0.079 


> 268 


66 


0.642 ± 0.079 



We will solve a direct problem and define the regions of Fj values which are compatible with equations (couplings) 

^i?,,xF,=J„ (i = 1,2,...,4) ^^^^ 
i 

where Ji is the observed flux of events with muon multiplicity from z-th interval - n^^ <i < n^(i_|_i) . 

Next we pass to the energy per nucleus and decrease the number of independent variables with the help of relations 

ME) - U{E) = U{E), (11) 

or 

h{E) = l.b!^{E), U{E) = f5{E). (12) 
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where fj{E) is the fraction of nuclei of kind j at the same energy per a nucleus. 

The relations (ITT]) are fulfilled at low energies {E,^ ~ 100 GeV) and the relations ([T^ are valid at i?„ — 100 TeV 
(see mass composition ([2])). We will find the solution of equations ([T0| in both cases, and in Section 7 discuss which 
variant f ([TT|) or is more preferable. 

Thus we work in the approximation fsiE) ~ fi{E) ~ fb{E) = 0.1467 and decrease the number of independent 
variables to three: /i, /2, /a- 

Passing from five variables Fj to three variables fk {k = 1,2,3), it is convenient to rewrite the equations pUj) as 
follows (we multiply and divide the j-th term by Aj-^ in each equation): 

3 

i?3y X Bj = J„ {i = 1, 2, 4) (13) 



where 



= Rrj/Ay, B, F, X Ay, J = 1, 2 

5 5 
j=3 j=3 



(14) 



In addition 



fk = BkX 



(15) 



To determine Bj we will use independent pairs of equations (|13p (for example for i = 1,2 or i = 2,3 etc.), and for 
closure of the equation system we use the normalization condition 

/1+./2 + 3/3 = 1, 

which (taking into account ([13]), (jlSp ) can be read so 

5 

B1 + B2 + 3B3 = E ^ 

As it is known, an inverse problem is incorrect (in our case the solution of system (jl3p is unstable at small variations 
of fluxes Ji). Therefore we solve the direct problem, namely: we define the regions of variables Bj which result in 
observed fluxes Ji to within one standard deviation ai = \/Ji. In the process we use the mass composition ([2|) as 
initial conditions. 

IV. PERMISSIBLE DOMAINS (I) 

In this Section we illustrate the procedure of determination of quantities Bj for the first two intervals of muon 
multiplicity: An^i = 114 - 151 and Anp2 = 151 - 189 (see table HIT)) . 
Let us write the equations (jl3l) for i = 1, 2 in an explicit form 

3 

^ ' -R3 1 j X Bj = Ji, 

(17) 

R32j X Bj = J2, 



or in a matrix form 



i?3_l X B = J_l, 



(18) 
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where 3x3 matrix R3_l is composed of elements R3ij for i = 1,2, and the third row of the matrix is the normalization 
conditions (fT6|) 

0.17374 0.38863 2.8885 
i?3_l = I 0.073225 0.16861 1.3701 | , (19) 
1 1 3 

the index _1 means that matrix R3 and vector J correspond to the first pair of equations (jl3p . Vector J_l is by 
definition (see table IIIip 

2.419 

J.l = ( 0.920 I (20) 
3.5453 

5 

(the third component of J_l is equal to the sum (ITB)) J.la = ^ Fj x A^'"^). 
With the help of ([131), dH]) we get 

0.88633 ' 

F = I 1.0991 I (21) 



0.51997 



where B™ are the initial conditions. 
Next we see that 

/ 2.0831 

i?3.1 X B*" = 0.9626 | (22) 
\ 3.5453 

To determine the regions of Bj, which satisfy the relations (|17|) we shall make the mass composition heavier (or 
lighter) B""'' until the result R x B™'' = J™'' faU outside the limits Ji ± cti, or J2 ± (72, where cti = \/Ji, 

o'2 = '\/i72 are the errors of a flux measurement. 

The change of the mass composition we shall realize by means of a decrease (or an increase) of the proton fraction 
Bi — > Bi =F (5 and an increase (or a decrease) of fractions of heavier nuclei B2 — ?> B2 ± (5/4, B3 — > B^ztS/A. Performing 
this procedure with a small step (for example S — 0.005 ) we define the regions 

B'""" < Bj < B^, 

which provide the fulfillment of the conditions (jl7p within one standard deviation. 
The conditions 

Jl - CTi < J™'' < Jl +(71, (23) 

and 

J2-(T2 < < J2+(T2. (24) 

define different regions of Bj , of course. As we work in the approximation of a slow change of the mass composition, 
one should choose an intersection of the regions as the solution for the mass composition averaged over the energy 
region under discussion ( 1170 < Epe < 2090 TeV, 2610 < Ep < 5840 TeV, see Fig. ^. 

Note the choice of the common region for solution of equations (jl7p means the use of two experimental points 
simultaneously. This reduces an influence of experimental errors. 

It should be noted that the regions of Bj values defined by may be disjoint at all. This would mean 

that either i) the mass composition is changed very rapidly (and our approximation Fj ~ const is incorrect), or 
ii) experimental data have such large errors which do not allow the simultaneous fulfillment of the conditions (|23p . 
(|24p . The latter variant is more plausible and this should keep in mind in what follows. 

Solving equations ([T7)l {i = 1,2) for the conditions ((^5)) we obtain (we present the results for variables fi, according 

to ^y. 

0.0808 < /i < 0.209 

0.274 < /a < 0.306 (25) 
0.172 < /a < 0.205 
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Figure 3: Energy distributions of protons and iron nuclei making a contribution to the flux of muon events with 114 < < 189. 
Areas under curves (p and Fe ) are equal to 1. The widths of distributions at half-height and fractions of events in these regions 
are indicated. 
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Figure 4: Permissible domains for the protons and He nuclei fractions according to the first pair of equations H13|) . i — 1,2. 



and for conditions 

0.207 < /i < 0.370 
0.280 < /2 < 0.321 
0.117 < /3 < 0.158 

The results are pictorially represented in FigH) The intersection of the regions (|25l) (pS)) is: 

0.207 < /i < 0.209 
0.280 < /2 < 0.306 
0.172 < /3 < 0.158 



(26) 



(27) 



As is obvious, the /s domains disjoint. We discuss a possible reason of that in Section IVIl As the temporary 
solution, we choose the average value - /s = 0.165. 

In a similar way, using the second pair of equations for i = 2, 3 (Anp2 = 151 — 189 and An^a — 189—268), we get: 
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for the condition J2 — cr2 < Jl"*" < J2 + (^2 



0.207 < /i < 0.370 
0.280 < /2 < 0.321 
0.117 < h < 0.158 

that coincides with (1^ , certainly, and for the condition J3 — (73 < J3 < 0/^3 + 0-3 

0.0808 < /i < 0.271 
0.278 < /2 < 0.326 
0.150 < h < 0.198 



(28) 



(29) 



The intersection of the regions (|28l) and (|29p (FigjS]) defines the mass composition in the range 151 < < 268 
(1630 < Epe < 3060 TeV, 3870 < Ep < 8750 TeV): 



0.207 < /i < 0.271 
0.280 < /2 < 0.321 
0.150 < /3 < 0.157. 



(30) 



0,5- 

He 

0,4- 



151 < n <268 



— ' 1 ' 1 ' 1 ' 1 — 

0,0 0,1 0,2 0,3 0,4 



Figure 5: Permissible domains for the protons and He nuclei fractions according to the second pair of equations (|13|) . i — 2,3 



Finally, the third independent pair of equations {i = 3,4, Art^3 = 189 — 268 and Anp4 > 268) gives the result ([29 
for J3 - CT3 < ^r'' < -^3 + 0-3, and for J4 - (74 < Jl"'' < J4 + 0-4 



0.192 < /i < 0.376 
0.279 < /2 < 0.325 
0.115 < /3 < 0.161. 

The intersection of the regions p9| and ((3T|) defines permissible domains for /, 

0.192 < /i < 0.271 
0.279 < /2 < 0.325 
0.150 < /3 < 0.161. 



(31) 



(32) 



The widths of energy distributions (at half-height) of iron nuclei and protons making a contribution to the flux of 
events with 189 < < 380 arc 2140 < Epe < 4270 TeV and 5330 < Ep < 12400 TeV. 



Note we used the integral point > 268 as Anf,4. In this case, R^j is determined by expression (|8}: Rij = Ijirif^ > 268). This has 



no affects on the final result. 
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V. PERMISSIBLE DOMAINS (II) 

Now we repeat the calculations of the previous Section using relationships faiE) — 1.5 x f4{E), J^iE) — fb[E) 
instead of ([TT|). 

Note that in this case 

i?3,3 = Ra/Al^ + '^(R^i/AY + iJ^sMs''). ^3 - y ^ x A]-\ (33) 

instead of and the normalization condition reads as follows (compare with p^ ) 

7 ^ 

Bi + + -5,3 = Y.F,x AY (34) 

(We assume h{E) = 1.5 x f^{E) i.e. Fa^P = 1.5^4^1-^, therefore factors 3/7 and 7/3 appear in ^ and ([M])) 
The matrix R3_l and vector _B™ take the form: 

0.17374 0.38863 2.1756 \ / 0.88633 \ 

ii3.1 = I 0.073225 0.16861 1.0268 , B"' = \ 1.0991 (35) 
1 1 2.3333 ) \ 0.66855 / 

Performing the procedure described it the previous Section, we obtain for j = 1, 2 

Jl - (71 < J™'' < Jl + CTl J2-(T2< Jl"'' <J2+(T2 

0.0554 < /i < 0.185 0.182 < /i < 0.349 

0.277 < /2 < 0.307 0.287 < /2 < 0.326 (36) 

0.230 < /a < 0.273 0.156 < /a < 0.211. 



The intersection of the regions (p6| is 



0.182 < /i < 0.185 

0.287 < /2 < 0.307 (37) 
0.230 < /a < 0.211. 



As with the approximation (fTTj) . the /a domains disjoint. We choose the average value /a = 0.220 as the temporary 

2 

3" 



solution. Note within the framework of our approximation fi = fz = I/a = 0.147. 



Permissible domains for i = 2, 3 i = 3, 4 are: 

i = 2,3 i = 3,4 

0.182 < /i < 0.250 0.160 < A < 0.250 

0.287 < /2 < 0.321 0.287 < A < 0.321 (38) 

0.203 < /a < 0.211 0.203 < /a < 0.218. 

VI. FITTING PROCEDURE 

Thus for the variant ^ f^{E) = 1.5 x fi{E), fi{E) ^ fz{E) (we shall call it "Model I" in the discussion that 
follows), we have obtained the following permissible domains of fi: 

1 = 1,2 i = 2,3 i = 3,4 



(39) 



0.182 < /i < 0.185 0.182 < /i < 0.250 0.160 < /i < 0.250 

0.287 < /2 < 0.307 0.287 < A < 0.321 0.287 < A < 0.321 

/a = 0.22 0.203 < /a < 0.211 0.203 < /a < 0.218 

fi = h = 0.147 0.135 < /4 = /s < 0.141 0.135 <fi^h< 0.145 

and for the variant ^ Jz{E) = fi{E) = /5(£;) ("Model 11") 

i = l,2 i = 2,3 1 = 3,4 

0.207 < /i < 0.209 0.207 < /i < 0.271 0.192 < /i < 0.271 

0.280 < /2 < 0.306 0.280 < A < 0.321 0.279 < A < 0.325 (40) 

/a = 0.165 0.150 <h< 0.157 0.150 < h< 0.161. 
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The obtained results have the status of permissible intervals of fractions fi. 

The second stage of data processing consists in narrowing of permissible intervals of fi. With this end in view, we 
carry out the simultaneous fit of 4 integral points (see table IIII[) and 3 finite diflference points: 

In, = I{> n^) , = 114, 151, 189, 268 , 

(41) 

J„^,_„;^ =/(>n^)-/(>n;,) , - n;, 114 - 151, 151 - 189, 189- 268. 

It will be now recalled that for i = 1,2 the domains do not have common region in both variants. This may 
be associated with experimental errors (see the remark before (|25p ) at a point (or at some points). It is clear from 
the general reasoning that the further from reality are experimental points, the smaller is the domain overlap. The 
smoothing procedure is needed in such a situation. 

The analysis has shown that the value of /(> 151) = 2.468 is incompatible with the assumption on a slow change 
of the mass composition in the energy range under consideration. In the present study we restrict the discussion to 
the correction of second integral point (/(> 151) — 2.468 ± 0.150) and the more exhaustive analysis will be presented 
in the later paper. Namely, we increase the second point value by ^cr: 

/(> 151) ^ 2.468 + 0.075 = 2.543. (42) 

This shift results in the common region of /a domains for i = 1 , 2 and in so doing the permissible intervals of fi 
become very close to the ones for z = 2, 3 and « = 3, 4. 
Next we perform fitting of 7 points (|4T|) requiring that: 

i) all data (7 points) are satisfied within the limits of experimental errors (±lcr), 

ii) fitted parameters (fi) are within permissible intervals. 

Taking into account the permissible intervals of fi (see (|39l) . (001)) take the form: 

Model I : 
0.182 < fp< 0.250 , 

(43) 

0.287 < fHe < 0.321 , 
0.203 < /jv < 0.211 , 0.135 < fst = fpe < 0.141 . 

Model II : 
0.207 <fp< 0.271 , 

(44) 

0.280 < fHe < 0.321 , 

0.150 < /at = fs, = fpe < 0.157 . 

In the process of fitting, fne and /^v have been chosen as independent fitted parameters. The rest fi are calculated 
with the help of expressions (|12p or (ITT|) . 

All possible sets of fi for the Model I are presented in Table IIVI For each admissible value of /jv was found the 
permissible interval of fne- In Table ITVl the mean and boundary values of fne are shown for each admissible /at. 
Columns Ij , Ji show the values of the fluxes (|41l) calculated at given values of fi {Ij and Ji have to be between min 
and max values which are shown in the second row) . All sets of fi presented in Table are statistically equivalent and 
can be written as follows: 

/„ = 0.2365 ± 0.0025 , fne = 0.2895 ± 0.0025 , /jv = 0.2038 ± 0.0008 , 

(45) 

fsz ^ fpe = 0.1356 ±0.0007 . 



^ The smoothing procedure gives the same results. 
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Table IV: Possible sets of primary nuclei fractions fi for the Model I 





fp 


fue. 


In 


fSi 


fpe 


{In A) 


/ll4 


/l51 


IlS9 


^268 


^114-151 


J^151-189 


Jl89-268 


min 


0.182 


0.287 


0.203 


0.1353 


0.1353 




4.678 


2.393 


1.427 


0.563 


2.199 


0.906 


0.815 


max 


0.250 


0.321 


0.211 


0.1467 


0.1467 




5.096 


2.693 


1.669 


0.721 


2.489 


1.084 


0.998 


0.333 


0.239 


0.287 


0.203 


0.135 


0.135 


1.929 


4.788 


2.581 


1.562 


0.707 


2.208 


1.018 


0.855 


0.313 


0.237 


0.290 


0.203 


0.135 


0.135 


1.933 


4.833 


2.605 


1.577 


0.714 


2.228 


1.028 


0.863 


0.333 


0.234 


0.292 


0.203 


0.135 


0.135 


1.936 


4.878 


2.629 


1.592 


0.721 


2.249 


1.038 


0.871 


0.317 


0.237 


0.287 


0.204 


0.136 


0.136 


1.937 


4.845 


2.612 


1.581 


0.716 


2.233 


1.031 


0.865 


0.359 


0.234 


0.290 


0.204 


0.136 


0.136 


1.941 


4.900 


2.641 


1.599 


0.724 


2.259 


1.042 


0.875 


0.339 


0.236 


0.287 


0.2046 


0.1364 


0.1364 


1.941 


4.879 


2.630 


1.593 


0.721 


2.249 


1.038 


0.872 



The fractions of light and heavy nuclei are 

fught =fp + fHe = 0.526 ± 0.005 , heavy = In + fs. + fpe = 0.474 ± 0.003 . (46) 

In the same manner, we obtain for the Model II: 

fp = 0.2405 ± 0.0045 , fne = 0.2985 ± 0.0145 , Jn = fs^ = fpe = 0.1535 ± 0.0035 . (47) 

hght ^fp + fHe = 0.539 ± 0.019 , fheavy = /jV + fs. + fpe - 0.461 ± 0.010 . (48) 

VII. DISCUSSION 

The procedure used at the first stage is based on the operation with two equations in 2 variables. To do this, it is 
necessary to set (in addition to a normalization condition) two relations between fractions fi. We use the relations 
([TT|) or (fT2|) . These additional relationships have a phcnomenological nature of course. 

Note the results obtained at the first stage depend on the initial conditions, but the second stage of data processing 
cancels this dependence practically. 

Let us remark also the method can be used under any energies (for example i?„ > 10^^ eV, > 1000) if the energy 
spectra of primary nuclei will be known. 

In regard to errors of fi determination, it should be noted the following. 

Experimental data (4 points of the integral spectrum and 3 points of the finite-difference spectrum constructed 
from them) are well described by a power function: 

I{n > n^,) = A ■ (114/n,J", J(n^) = A ■ (m/114) • (114/n^)(™+i) 

(49) 

A = 4.9442 ± 0.1357, m = 2.3695 ± 0.0921. 

As we can see, the relative error in spectrum parameters is about 4%, while the relative error of the initial experimental 
data on the integral spectrum is about 5-10%. This decrease of the relative error is due to matching data with the 
function and its derivative to get the results p9)) (A and m), while the initial errors of 5-10% are applicable only to 
the function. 

Taking into account (|49p . we are sure that correct processing of the experimental data must detect the mass 
composition with relative errors in fractions fi under 5%. Of course, these errors can be bigger due to additional 
simplifying assumptions, e.g. MODEL I or MODEL II. In any case, if errors over 5% are indicated in the final result, 
then the precision of the method used to process the experimental data (i.e. the solution algorithm for the ill-posed 
problem) is lower than the precision of these data. 

We now discuss the results of the paper. Numbers after ± signs in (H5t and (|T7)) are not the fitting errors. These 
formulas are compact descriptions of restricted domains, whose full description is given in Tables. As for the true 
errors of fi, the situation is as follows. The errors can be computed only after determining ALL solutions, matched 
with the experimental data via integral and finite difference spectra. The used algorithm is mathematically very rigid, 
since it requires matching independent solutions of three pairs of equations with the solution of seven equations. This 
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algorithm does not detect all possible solutions but only a part of them. But the algorithm itself is compatible, so it 
is possible to claim that the permissible intervals of fi are expanded by at most 5% of the determined ones. If we also 
take into account ambiguity of the interaction model choice (appendix A), then we get the error value about 10%. 
In this context MODEL I and MODEL II result in the same results: 



fp = 0.236 ± 0.020, fHe ^ 0.290 ± 0.020, fn = 0.474 ± 0.030 (50) 

In our opinion, MODEL I is more preferred for the following reasons. First, the relations ()12|) are satisfied for mass 
composition © defined a.t E ~ 100 TeV, i.e. in the neighboring energy region. In the second place, the same relations 
are characteristic for chemical composition of shells of II type supernova stars, which (according to current concepts) 
are the sources of high energy cosmic rays. 

The results presented above are shown in FigEl One can see that uncertainty intervals of fi values for the Model 
I are noticeably less than the ones for the Model II. Possibly this is one more circumstance pointing to the most 
adequacy of the Model I. 
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Figure 6: Permissible domains for primary nuclei fractions fi. 



VIII. CONCLUSION 



We have presented the method of retrieval of information on CR mass composition based on a solution of the 
direct problem ()10|) with a determination of permissible intervals of primary nuclei fractions fi. At the second stage, 
we constrict the intervals of fi with the help of fitting procedure using the information about the integral spectrum 
I{> n^j) and its derivative J{An^). 

Data processing is based on the theoretical model representing the integral muon multiplicity spectrum as the 
superposition of the spectra corresponding to different kinds of primary nuclei. In so doing, it should be kept in mind 
that we have fixed the interaction model (see appendix A). 

With the framework of three components (protons, helium and heavy nuclei) and under the assumption on a slowly 
change, the CR mass composition in the region 10^^ — 10^^ eV has been defined: 

fp = 0.236 ± 0.020, fHe = 0.290 ± 0.020, fn = 0.474 ± 0.030 (51) 

The result ((5T|) should be read as the estimation (rather precise) of CR mass composition in the energy region of 
10""^^ — 10^^ eV. Thus our analysis points out that CR mass composition become some more heavy in comparison with 
the one 

The method can be used under any energies if the energy spectra of primary nuclei are known. 
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IX. APPENDIX A 

To calculate the parameters A{m) and G(m) [H, we have used Monte-Carlo shnulation and the approximation 
of spatial-energy distribution function (SDF) of muons in EAS obtained in (22j . 

/(r, >E,Eo)=C^ exp[-{r/ro)\ {A.l) 

here 

0.95 0.42 , 0.2 



° (1 + i2.5£:)o-92 £:i.23£;o.9' ' 0.2 + £;„ ' 

Eo is an energy per nucleon in a primary nucleus, E is the muon threshold energy (E in TeV, r in meters), C - a 
normalization factor. 

The muon density p at a distance r from EAS axis is difined by the expression 

p(r, > E, Eo)2TTrdr = ri^{A, E^, > E)f{r, Eo, > E)rdr, {A.2) 

here is the average number of muons with energy > E produced by a primary nucleus with energy _B„ = AEq (A 
is the number of nucleons in a nucleus) [2^ : 

0.0187^(6*)^ /£;o\"-^V Eo 



^M,Eo,>E,0) = —^^-f) (A3) 
where Eo and E in TeV, 

r..,,r.^ f r. 1 + 0.36 X /n(cOS 6*) 

lg(W + U.5Eo) cosd 

- zenith angle. 

X. APPENDIX B 

The negative binomial distribution is a discrete probability distribution of the number of successes in a sequence 
of Bernoulli trials before a specified (non-random) number k of failures occurs 

B(fc,p) = C7;+'=-i(l-p)V, fc = 0,1,2,... (B.l) 

here p is the probability of success. 

Putting r = and taking into account = kp/{l — p) we obtain 

" \n^,/k + \J (n^/fc + 1)'^ 
The parameter k was chosen in the form obtained in [2^ 

X 0.18 



/ = A 1 -t- 0.013-^ (B.3) 



p-r ' V E 

At such choice of k the variance of B(n^,n^, k) distribution is greater than the one of Poisson distribution in 6 — 10 
times for protons and 2 — 3 times for iron nuclei in the energy region under discussion in accordance with results of 
others works (e.g. [13, Hi]). 
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